#####################################################################################################
# This script generates a map of disenfranchisement in antebellum Kentucky                          #
# It draws on existing data series on the location of lakes, rivers, canals, and railroad,          #
# as well as the topography of the state.                                                           #
#                                                                                                   #
# Topography data is shown, but not necessary for replication. This can be commented out if you     #
# have difficulties with get_elev_raster.                                                           # 
#                                                                                                   #
#####################################################################################################

### Start by loading packages used across different maps. 
rm(list = ls())

library(shapefiles)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(ggplot2)
library(readxl)
library(dismo)
library(raster)
library(rgeos)
library(elevatr)
library(ggmap) 
library(tibble)

### Provide version list. 
sessionInfo()

### Set the API key to get the Google geocoded location of cities. For users without a key, longitude and latitude are provided below
#register_google('[API key here]')
### Set the directory
dir <-paste("../", sep="")
setwd(dir)

### Create a heat map function for maps. 
plot.heat <- function(states.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(states.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(states.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  states.map@data$zCat <- cut(states.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(states.map@data$zCat)
  if (is.null(col.vec)) col.vec <- color.palette1(length(levels(states.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(states.map@data$zCat) <- cutpointsColors
  plot(states.map,border=gray(.8), lwd=bw,axes = FALSE, las = 1,col=as.character(states.map@data$zCat))
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
}

### Define color scheme and create index for figures
color.palette1 =  colorRampPalette(c("black","white"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette1(6)
index = c(cols[1], cols[3], cols[4], cols[5], cols[6])

################################################################################
###                        Kentucky Census of Voters                         ###
################################################################################

### Set figure details and year range
year <- paste("1850")
year2<- paste("1860")
lgb1		= "60-70%"
lgb2		= "70-80%"
lgb3		= "80-90%"
lgb4		= "90-95%"
lgb5		= "95-100%"
lg1 		= "650"
lg2		= "300"
lg3		= "100"
lg4		= "50"
lg5		= "1"
label <-paste("Figures/Supplemental/Supplemental-Figure-3.pdf", sep="")

## Create list of inputs for shape/voter files. One shapefile, for Kentucky in 1850. 
## One voter files

input1 <- paste("Shapefiles/Kentucky", year, ".shp", sep="")
input2 <- paste("Data/KY-Census-Of-Voters.csv", sep="")

## Load rail, rivers, canals, elevation, and the county level data. 
rail <- readOGR("Shapefiles/RR1826-1911Modified0509161", "RR1826-1911Modified050916")
rail <- rail[rail$InOpBy < year,]
rivers <- readOGR("Shapefiles/SteamboatNavigatedRiverContigUSAprojection", layer = "SteamboatNavigatedRiverContigUSAprojection")
rivers <- rivers[rivers$START < year & rivers$END > year2,]
canals <- readOGR("Shapefiles/19thC_Canals_March2017", layer = "19thC_Canals_March2017")
canals <- canals[which(canals$OPENED < year & canals$CLOSED > year),]

counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
counties.map <- spTransform(counties.map, CRS(proj4string(rivers)))
canals <- spTransform(canals, CRS(proj4string(rivers)))
rail <- spTransform(rail, CRS(proj4string(rivers)))
elevation <- get_elev_raster(counties.map, z = 8, expand = 150000)

data.counties <- read.csv(input2)

### The location of the cities requires user to have a Google API key. 
### I comment out the command accessing Google's geocoding, and replace it with the longitute
### and latitude of the relevant cities.
cities1 <- c("Louisville, KY", "Covington, KY", "Frankfort, KY", "Lexington, KY", "Maysville, KY", "Paducah, KY")
#cities <- geocode(cities1)
cities <-tibble(lon = c(-85.7585,-84.5086,-84.8733,-84.5037,-83.7444,-88.6000), lat= c(38.2527, 39.0837, 38.2009, 38.0406, 38.6412, 37.0834))
### Input estimated disenfranchised population of cities, from county information and census records.
cities2 <- c(4350, 705, 342, 244, 164, 152)
### Remove state abbreviation for map
cities1 <- c("Louisville", "Covington and Newport", "Frankfort", "Lexington", "Maysville", "Paducah")
cities <- SpatialPointsDataFrame(cities, data = cities, proj4string = CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
cities <- spTransform(cities, CRS(proj4string(rivers)))

names(data.counties )[2]<-"FIPS"
names(data.counties )[5]<-"est"

counties.map@data$order <- 1:nrow(counties.map@data)
counties.map@data <- merge(counties.map@data,data.counties, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.map@data<-counties.map@data[order(counties.map@data$order),]
names(counties.map@data)[ names(counties.map@data)=="est"]<-"noise"

### Create Supplemental-Figure-3 
pdf(label, width=8, height=6, bg="transparent", family="sans" )
par(mfrow=c(1,1),oma=c(0,0,0,0), mar=c(0,0,0,0))

plot(counties.map,density=c(30))
par(new=TRUE)
plot.heat(counties.map,z="noise",breaks=(c(60, 70, 80, 90, 95, 120)), plot.legend=FALSE)
plot(rail, col = "indianred", add = TRUE)
plot(rivers[rivers$END > 1860,], col = "dodgerblue", lwd = 1, add = TRUE)
plot(canals, col = "blue", lwd = 1, add = TRUE)

### Create the topographical shading. This is not essential to the replication. 
### If user has difficulty getting elevation data from get_elev_raster above, then
### section can be commented out.
mosaic <- elevation * 10 
slope <- terrain(mosaic, opt="slope", unit='radians')
aspect <- terrain(mosaic, opt="aspect", unit='radians')
hillshade <- hillShade(slope, aspect, angle=45, direction=315)
hillshade2 <- aggregate(hillshade , fact = 5 , method = "bilinear" )
hillshade2 <- focal(hillshade2, w=matrix(1/9, nc=3, nr=3), mean)
slope2 <- aggregate(slope , fact = 5 , method = "bilinear" )
slope2 <- focal(slope2, w=matrix(1/9, nc=3, nr=3), mean)
plot(hillshade2, col=grey.colors(100, start=0, end=1), legend=F, alpha=.1, add=TRUE)
plot(slope2, col=grey.colors(100, start=1, end=0), legend=F, alpha=.1, add=TRUE)

### Plot location of cities / disenfranchised adult free men.
plot(cities, add = T,  pch=21,  bg=alpha(rgb(1,1,1), 0),col="black", lwd=.4, cex=cities2/300)
text(cities, labels=cities1, pos=3)

### Insert index
legend('topleft', bty="n", fill=index , bg="transparent", ncol=1,
       x.intersp = 1, xjust=0, yjust=0, text.width=10, cex=.9,
       legend = c(paste(lgb1, "Adult White Men\nEnfranchised"), lgb2, lgb3, lgb4, lgb5))

dev.off()

